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ABSTRACT 


A computer simulation model of the phosphate, phytoplank- 
ton and zooplankton dynamics in Monterey Bay was examined and 
modified. The model is driven by four foreing functions ex- 
pressed as annual cycles of upwelling velocity, incident 
Solar radiation, mixed layer depth, and mixed layer tempera- 
ture. An alternate upwelling index was developed based on 
the local wind field. A revised radiation index is employed 
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Analysis of the model's response to Sinking and advection of 
phytoplankton was examined. The importance of seasonal in- 
creases in predators was introduced as a controlling factor 
in the seasonal growth of zooplankton. The model is able to 
predict the seasonal trends of phosphate, phytoplankton, and 
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Lo EGRODUCTRON 


Pe CURE OSE, 

In the interests of assessing the biomass of oceanic 
areas, various mathematical models have been designed to 
predict the dynamic response of ecosystems to change in the 
physical and chemical environment. Such models necessarily 
rely on accurate characterization of the forcing functions 
and accurate representation of the system with equations 
which are consistent with the specific region under study. 

Recent attention paid to modeling systems in upwelling 
regions, Coastal Upwelling Ecosystems Analysis (CUEA),has 
Drought about developments in both these areas. Adequate 
data of a physical and chemical nature is becoming available 
and significant progress in the refinement of governing equa- 
tions 1s evident (Walsh, 1973). 

The research presented here was aimed at creating a simu- 
lation model to describe the seasonal plankton dynamics in 
Monterey Bay, California, as known from the best available 
data. An existing simulation model (Pearson, 1975) was 
evaluated in an effort to make refinements and to judge el 
applicability of the time simulation technique to the 


Monterey region, characterized by seasonal upwelling. 


B. @SM@DELING PHILOSOPHY IN°THE ECOLOGICAL CONTEXT 
The merit of mathematically modeling systems of a purely 


physical nature in which parameters, initial conditions, 
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boundary conditions and the variable relationships are readi- 
ly specified 1s unique. The value of current models of a 
biological nature requires some comment. 

According to Odum (1971) there are a number of reasons 
behind constructing any mathematical model. Prediction of 
future states of the variables involved is an end in itself. 
Use of the prediction accuracy to focus data acquisition 
efforts or direct attention to deficiencies in concepts 1s 
likewise valid. If the model is real in the sense that it 
attempts to predict changes in the state variables based on 
a framework drawn from research in the real world, then 
failures may be useful in delineating weaknesses in the 
governing equations. 

Construction of a "real" model is at best difficult. A 
perfect model of an ecosystem assumes that "all states and 
rates of change of the variables are known at all times" 
(Walsh, 1973). Given the vastness of complex interactions 
in even the simplest of biological systems, it may be stated 
that "no perfect representation of the real world exists 
emeept the real world itself" (Walsh, 1973). 

It follows then that certain assumptions, logical state- 
ments and approximations of real dynamics are necessary to 
allow the model to operate. The resulting model may bear 
little correlation to the actual physical, chemical and 
biological situation at hand. Again, the argument for con- 
structing the model must be considered. For purposes of 


prediction, a model which makes a logical statement relating 


it 





two state variables in the attempt to specify a relationship 
which may be unknown in part or in whole has merit. The 
Tact that the logic of the model is a crude simplification 
is incidental if the prediction capability is proven. 

All models can be characterized on the basis of realism, 
precision, and generality (Odum, 1971). It has been suggested 
that these three qualities are mutually exclusive in ecologi- 
cal models (Patten, 1971). A fourth category, simplicity, 
might be considered. Exacting detail, while lending to 
realism and precision often forces the model to be specific. 
Only when the model becomes an isomorph, 1.e., a one-to-one 
correspondence to the real world, will generality be restored 
(Weish, 1973). 

The problems in creating a realistic model have been dis- 
cussed briefly. The value of this type of model lies in re- 
search guidance. Where a specific question is asked of an 
ecological model, precision in prediction capability may be 
enhanced by sacrificing the other qualities. Such an 
approach might be taken in a fisheries model where a precise 
Output for a limited region is desired (Odum, 1971). Simplic- 
ity and generality at the expense of both precision and 
realism may enhance understanding of broad ecological con- 
@eDts. 

The study presented herein was conducted with the idea 
of creating a model having good predictive capabilities for 
conditions occuring over a long span of time in a limited 


oceanic region. 
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Ceo wee eR DP TTON OF THE MODEL 

The original ecosystem model was developed by Pearson 
(1975) as part of continuing research into the correlation 
or zooplankton biomass with the chemical and acoustic pro- 
perties of the ocean (Traganza, 1976). It is the intent of 
this research to identify areas of weakness in the model 
and make the appropriate changes to the equations and input 
parameters. These changes include a refinement of the up- 
welling index and radiation index as well as inclusion of 
advection, sinking and predation terms in the governing 
equations. 

The modeling technique uSed is a time dependent simula- 
tion solved digitally on an IBM 360 computer using the IBM 
Continuous Simulation Modeling Program (CSMP-360). The 
dynamic equations are written in non-linear differential 
form and driven by four exogenous forcing functions expressed 
as annual cycles of upwelling velocity, mixed layer depth, 
mixed layer temperature, and incident solar radiation. The 
output consists of annual cycles of phosphate concentration 
(microgram-atoms of phosphorous per liter, ug-atP/2), phyto- 
plankton biomass (grams of carbon per square meter, gC/m*) 
ane Standing stock of herbivorous zooplankton (g C/m*) ng 
the mixed layer. 

The Pearson version (Figure 1) of the model is defined 


by the following basic equations: 
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Pictoral representation of Monterey Upwelling 
Ecosysten™ 
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where: 
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UPTAK 
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dxX2 
at 


PROD 


RESP 2 


GRAZ 


dx3 
ie 


GRAZ 


RESP 3 
VOID 


LOSS 


time rate of change of phosphate concentra- 
tion (yug-at P/1) in the mixed layer 


input of phosphate from upwelling and mixing 
from below the mixed layer (ug-at P/1) 


phosphate recycled as biological excretory 
products within the mixed layer (yg-at P/1) 


phosphate depleted by phytoplankton utiliza-~ 
ttom (Ug-=at P/1) 

PROD -PRESP2 — GRAZ ©) 
time nate of change of phytoplankton biomass 
CeCym-= ) 

DNOtOsynunetLle™ proauer? on (gC/m*) 


respiration losses (gC/m*) 


‘grazing losses (gC/m*) 


GRAZ - RESP3 - VOID - LOSS Gap 
input, due to ingestion of phytoplankton 
Cacyim™ ) 


respiration losses (gC/m*) 
excretory losses (gC/m*) 
Beck depletion =) mori lity (gC/m*) 


(Pearson, 1975) 
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Il. BACKGROUND THEORY 


Be. UPWELLENG THEORY 
1. Ekman Dynamics 

Seasonal upwelling of nutrient-rich deep water into 
the productive mixed layer is an important aspect of the 
-Monterey Bay ecosystem. The principal source of phosphate 
in the mixed layer is vertical advection which is produced 
in response to wind-induced Ekman circulation. Theoretical 
calculations of the offshore component of horizontal current 
are based on Ekman pure drift theory in an infinitely deep, 
homogeneous ocean (Neumann and Pierson, 1966). The compo- 


Memes Of the horizontal current take the form: 


-(1/D)Z 
e 


Ce Ns COS[(45-(9/D)Z)] and (4a) 
vez vy e7$7/P)4 sin cus-(m/D)Z (4b) 
D = (A/pw SIN oe and (Sia) 
ve = t,/Aa Cab) 
where: 
U = velocity component in the x direction 


¥ = veloclty coOmpement in “the y direction 
(parallel to wind) 


V. = surface velocity component (45° to ease hiceeone 
wind, Northern Hemisphere) 


= Cepen 
D = Ekman depth of frictional resistance 


A = coefficient of eddy viscosity 
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tT = surface wind stress in y direction (dynes) 
a = fo/A 
@ = latitude 
op = density of water (g/em*) 


Ff = Coriolis parameter = 2u SIN ? 


The coefficient of eddy viscosity, A, was assumed 
constant Since detailed information on the small scale mo- 
tions of the surface layer was not available. This seats 
tion is generally made in mass transport studies (Neumann, 
196s ). 

The mass transport due to wind driven current is 
found by integrating the equations of velocity (U and V) 


Over the depth of the water column, as shown: 


cm == Ua = t./Ff (6a) 
x - y, 
and = = 9 VdZ = 0 (6b) 
where: 
Sy = Mees transport in the x directzon 
Sy = mass transport in the y direction 


From the equations for mass transport, it is observed that 
the net transport is 90 degrees to the right of the wind in 
the northern hemisphere and that it is directly proportional 
to the wind stress acting on the sea surface, parallel to 


the coast. 


Ly 





2. Wand Stress on the Sea Surface 
The coupling of wind energy to the water at the air- 
sea interface is defined by the conventionally-accepted form 


of the wind stress equation: 


yy = oC We (Gordon, 19172) . CH) 
where: / = wind stress parallel to the coast (dynes/em*) 

o' = density of air 

C, = drag coefficient at height z (non-dimensional) 

ie = wind speed at height z (cm/sec) 


As in other momentum exchange studies, the dynamic 
coefficients appear to present the singular most difficult 
problem. Wilson (1960) states that the drag coefficient for 
air flow over a water surface is dependent on the wind speed 
(W), the height of measurement(z), a surface roughness para- 
meter (20); and atmospheric stability terms. There may be 
additional dependence on oceanic parameters of depth, fetch 
distance and wave height. An average value of the drag co- 
efficient for high wind speeds (i.e. greater than 10 m/sec) 


a ioe se BO 


is arrived at by examination of research com- 
pleted through 1959 (Wilson, 1960). 
Work by Deacon (1962) established an empirical rela- 


tionship for winds up to 13 meters/second which gives a 


linear dependence of C_ (at z = 10m) to the wind speed as 
rollows: 
a 
Cig = Cie + Oo Wig? x 10 for W in m/sec (8) 


(Roll, 1965) 


18 





Deacon's research was carried out using data from ship obser- 
vations and coastal regions. 
3. Spatial Extent 

The offshore movement of water in a region bounded 
by a coastline causes replacement water to be upwelled from 
below the layer affected by the surface wind stress. Accord- 
ingly, it is necessary to specify the spatial dimensions of 
the region in order to apply principles of continuity in de- 
termining the vertical current speed. The significant para- 
meters involved in upwelling are: (1) the offshore horizontal 
dimension of the upwelling region, which is most directly 
related to the width of the wind field; (2) the offshore 
component of the Ekman current (U); and (3) the depth of the 
Ekman layer (D). 

Estimates of the maximum depth of the source of up- 
welled water in a coastal region are approximately 100 to 
200 meters (Neuman and Pierson, 1966) based on the slope of 
the isotherms. The mass transport offshore due to wind 
drift current has been previously detailed as occurring in 
a surface layer of depth, D, corresponding to the Ekman depth 
of frictional resistance; regardless of the source depth of 
the upwelled water, it transgresses the "boundary" at depth, 
D, before being carried offshore. The vertical dimension of 
the region is then readily specified. It is assumed that 
the depth of frictional resistance (D) and the depth of the 
mixed layer (Z) are approximately coincident. The differ- 


ences between (D) and (Z) are important when phytoplankton 
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are advected out of the Ekman layer but exist throughout 
=heemaxed layer. 

The horizontal extent of the upwelling region 1s 
considerably more nebulous than the depth of the mixed layer. 
Classical estimates of this dimension limit the zone to 
probably no more than 100 kilometers in offshore width. It 
may be expected that the width of the region is dependent 
on several factors; among them, the width of the wind field, 
variations and non-linearities in the energy exchange pro- 
cesses due to surface roughness characteristics or stability 
changes in both the atmosphere and the water column, as well 
* Spatial and temporal oscillations of the significant wind 
weetor. In addition, local coastal and bottom topography 
may figure extensively in the problem as will be discussed 
later. 

Sverdvrup (1938) observed a relatively well defined 
offshore boundary to a coastal upwelling zone, coincident 
Witha downwind current which was marked by an intense verti- 
@aml Gradient of velocity (Sverdrup et al., 1942). He further 
observed an offshore migration of the current band at a rate 
somewhat less than the speed of the offshore surface current 
within the upwelling region. The latter fact implied the 
possibility of cellular circulation patterns in the near sur- 
face upwelling zone. 

Hidaka (1954) proposed a steady state theory of up- 
welling in which he defined a horizontal frictional distance, 


D,, analogous to the depth of frictional resistance 
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appearing in Ekman's work. Applying reasonable values to 
Hidaka's expression confines the total cellular circulation 
region to a width of about 15 kilometers. The width of the 


region is defined by: 


L, Phe i 
Dy = 1(2A,/2w SIN ¢)* where: Ay = 10’ (gem Sees 


) cD 
(Smith, 1968) 
Downward vertical currents occurring in the seaward half of 
the cell limit the upwelling to a width of 7.5 kilometers. 
The theoretical width 1s not consistent with Eee eaee of 
upwelling at greater distances from the coast (Barnes, 1969). 
Yoshida (1967) developed a quasi-steady state model 
applicable to eastern boundary current regions (Smith, 1968). 
An expression was derived for the horizontal dimension of 


the coastal upwelling region which is given by: 


- 


2 
2 L 
L. = eae = (Hu/Dy)* for latitudes greater than 
$ 7 22 degrees (10) 
where: Ly = horizontal width of the coastal boundary 


region (m) 


= gravitational acceleration cpm +) 


& 
H = thickness of the upper layer (m) 
D = H + thickness of the lower layer (m) 
F = Coriolis parameter (sec +) 
0 = density of water (g-om™ 7) 
Ao = density difference between deep and surface 
layers (g-cm73) 
“uo = Intermal trictien @oett teres jee) 
vid eee xX ioe wece (Smith, 13m) 
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B. CURRENTLY AVAILABLE UPWELLING INDEX 

Bakun (1976) has calculated coastal upwelling indices 
for the coast of North America using Fleet Numerical Weather 
Central monthly mean pressure fields to compute the geos- 
trophic wind. Analysis is done on a 63 x 63 point grid of 
approximately 200 nautical mile mesh length (Bakun, 1973). 
Figure 2 shows Bakun's results for 1974 at 36 degrees North, 
122 degrees West, approximatley 50 miles south of Monterey 
Bay. In reviewing several aspects of hydrographic surveys 
taken in Monterey Bay, it was suspected that this region is 
characterized by local upwelling patterns not appearing in 
the index computed by Bakun. Further discussion follows in 


mae methods section. 


@ fFrrecl’s OF UPWELLING ON INCIDENT RADIATION 

Cold water upwelled from depth during the spring and 
summer months brings about both advection fog and stratus 
cloud formation. The fog occurs as warm surface air 1s 
cooled by the cold seawater to its saturation point causing 
the condensation of the contained water vapor. Stratus 
clouds occur below the base of a quasi-permanent atmospheric 
temperature inversion (Tont, 1975). 

The net effect of both fog and stratus layers 1s to re- 
duce the amount of solar radiation reaching the sea surface 
during the upwelling period. Tont (1975) obtained mean 
values of solar radiation at the surface over the period 
1950-1973 at San Diego. A correlation study of the annual 


upwelling index and the amount of sunlight at the surface 
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showed a maximum of about 40 per cent reduction in available 
sunlight during the peak upwelling months. The result of 


Tont's research is shown in Figure 3. 


D. SINKING OF ALGAL CELLS 

Review of sedimentation rates of algal cells indicate 
justification for including a sinking term in the phytoplank- 
ton equation. The Pearson (1975) model assumes uniform algal 
concentation within the mixed layer. It 1s reasonable that 
changes in the concentration of phytoplankton caused by ver- 
tical movement out of the layer will occur. 

Parsons and Takahashi (1973) state that the theoretical 


Sinking rate of phytoplankton can be determined from: 








0 
vy = 28h (Pi=8) (11) 
9 ae 


where: 

r = radius of the cell 

g = gravitational acceleration 

o0' = density of the organism 

o > density of the medium (sea water) 

fl =ovaccOcaty Olan medium 
6 = form resistance coefficient (accounts for non- 
spherical shapes) : 

Measured rates of sinking for live algal cells vary be- 

tween zero and 30 meters per day according to Parsons and 


Takahashi (1973). These rates represent a broad range of 


species and sizes. 


24 





aeironts SJtalosoud INasvdad 


"(GL6T “3UOL 
goe1je) euTysuns eTqtssod jo jUus0ued puP 
ic _-7a) xoput BSutttemdn jo sueveu ATYWUOW °E 2ANBTY 


i N39 Ad ---o -- 
ONVITS Mell —o=— 





25 





Be wADVECTION 

Dimi Knew thet salamafieamt current velocities occur 
within the — from Ekman dynamics. It might be ex- 
pected then that the commenymepaien of phosphate, phytoplank- 
ton and zooplankton will be affected by advection. O'Brien 
and Wrobleski (1972) have shown by scale analysis that advec- 
tion and biological productivity of the Peru upwelling eco- 
system are equally important. 

The advective terms of the material derivative of a 


property (C) are written: 


dC dC dC (12) 


Advective change = U ax V — * W aa 
where: U,V,W = velocity components in x, y, and z directions 
respectively 
ac oC | AC 


a%? 9% 97 7 gradients ef the concentration of property 
(C) in the x, y, and z directions 
A review of California Cooperative Fisheries Investiga- 
tiene (GALCOFI) data for 1974 (CALCOFI, 1974) shows that the 
phosphate concentration varies significantly in the vertical 
direction but the horizontal variation is not as significant. 
It is assumed for the purposes of this model that the ad- 


vection term for phosphate reduces to the vertical component: 


Advective change = W a7 nay 
oe axl , aXl _ 4 
sanee: ax ay 
Pearson (1975) has included this term in the computer model 
DXl - XL 
° = (i So ee 
as: UEREEL W =XL (14) 
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where: 


W = vertical speed (m/day) 


EKL = depth of the Ekman layer (m), EKL = D 
DX1 = concentration of phosphate below the mixed 
layer (ug-at P/1) 
Xl = comeenaration of phosphate in the mixed layer 
(constant) (yug-at P/1) 
a = vertical gradient of phosphate (ug-at P/m 1) 


The advection of zooplankton in the revised model is 
discounted because it is assumed that they have developed 
mechanisms which permit them to remain in a particular 
oceanic region by swimming or riding cellular circulation 
evmrents. 

Phytoplankton, however, are affected by advection. 
SemeeOrk data (CALCOFI,,.1974) shows that the significant 
gradients of phytoplankton occur in the vertical and off- 
shore directions; therefore, the advection expression for 


phytoplankton is reduced to: 


Advective change = U =~ + W == CES) 


The horizontal velocity is related to the vertical veloc- 


ity by continuity (when it 1s assumed that oad = 0, parallel 


to the coast) such that: 


U = - WL/D GuGy) 


2/7 





where: = horizontal velocity component offshore (m/day) 
upward vertical velocity component (m/day) 


= width of the region (m) 


O a =; = 
i 


= Ekman depth (m) 


The phytoplankton advection term can be written: 


Advective change = W [- = a i a] Gy 


The gradients of phytoplankton can be approximated by: 


G2. LIM AxX2 OxX2.. Jbl aye7 
ss and —=— 


= «(AXSO AN a” he he (18aéb) 


If it is assumed that the concentrations of phytoplankton 
below the mixed layer and outside the modeled area are much 
less than the concentration in the modeled region, then the 


gradient can be approximated by: 


- 


X2=-0 


= horizontal gradient; (19) 
at - vertical gradient; G20) 
where: X2 = phytoplankton concentration in the 


mixed layer 


AX = horizontal distance at which the concen- 
tration of phytoplankton becomes negli- 
gible with respect to the modeled area 


AZ vertical distance at which phytoplankton 


concentration becomes negligibly small 
The phytoplankton advection term can then be written as: 
Advective change to L 1 


Phytoplankton in the = W(X2) [- pay + 77! 
modeled area 


C2a@) 
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The sign convention is established as positive (+) into the 


model area and negative (-) out. 


e, PREDATEON PRESSURE 

Under constent environmental conditions, the dynamics of 
zooplankton growth are a function of the food supply and 
predation pressure. Natural processes of population control 
insure that a given species will not be eliminated by low 
food supply and high predation pressure or that unstable 
growth will not result from the inverse situation. Attempts 
to translate the stabilizing factors into mathematical terms 
often result in simplified general relationships that do not 
exhibit the flexibilities inherent in the real system. Not- 
Withstanding the limitations, a predation term is included 
in the simulation model. 

The long term survivability of both predator and prey is 
keyed oerly to maintenance of a balance in energy expendi- 
tures and gains. To date, laboratory experiments to dupli- 
cate this balance have not been entirely successful. Prey 
stocks in small laboratory environments have been artificial- 
ly supported by providing refuge where predator access to 
the prey was denied and by introducing additional prey to 
hie system to replenish the supply (Patten, 1971). Neverthe- 
less, population control by predation is prey-density depen- 
dent as a first approximation. For a system of one prey 
species and a single predator species, predation pressure 
may be represented by: 

Ken D tenets. D<d- and 


1 I 
fomeped 


'U 
iN 


k (22) 
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where: P = predation pressure or food intake per unit 
predator 


D prey density 
d,= satiation threshold of prey density 


Ky and Ko = constants 


Above prey denoaty vabues»of dj > predator satiation will oc- 
cur and pressure will level off at some Conereant Salah vale, 
The function is shown on an arbitrary scale in Figure 4 
according to Patten (1971). This simple relationship is 

not entirely satisfactory for the predator since low food 
Supplies would result in starvation. Under such circum- 
stances, and given a second prey species of perhaps less 
Palatability, the predator would switch food intake to the 
more abundant species of prey. There are several obvious 
implications to the switching phenomenon: 

1. predation on the most abundant species assures 

stable growth of that species, 

2. Switching from a declining species serves to prevent 

elimination @@ that spécres (Paittten, 1971), and 

3. elimination of the predator by overgrazing one avail- 

able food source is avoided. 

The model under consideration in this thesis is a single 
species approximation in which predation pressure is incor- 
porated with other zooplankton mortality factors, l.e. 
morphological death, in the "LOSS" term. This term is given 


by Pearson (1975) as: 


LOSS = Lie G33 (230 
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Figure 4. Predation pressure as a function of prey 


density (zooplankton biomass) (after 
Patten, 1971). 
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where: LOSS = fraction of the zooplankton population re- 
moved from the system in a day's time 
(g C/m2 day) 


b =S6onstane rate of less of percenr, loss pen 
day (day71) 


o-oo p Lani tones! Stand tne crop. a Ge C/m?) 


Harriston concluded in 1960 that "herbivores are preda- 
toe immited™ (Patten, 1971). 

Conclusions by Patten (1971) show that herbivores are 
both food and predator limited. Therefore, in the modified 
model, the LOSS term is reserved for variable predator limi- 
tation of herbivorous zooplankton and natural death of 


zooplankton herbivores is assumed to be negligible. 
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itl. 2eieoDs 


The existing model appears weak in several areas. Orig- 
inally, the forcing functions of upwelling and solar radia- 
tion which are critical to the simulation were obtained from 
average cycles for the California coast (upwelling) and the 
latitude band between 30 and 40 degrees north. Phytoplankton 
Sinking and advection terms were lacking and the predation 


pressure term in the zooplankton equation needed improvement. 


A. UPWELLING PROGRAM 

A computer program was written to calculate the annual 
upwelling index for Monterey Bay. Wind data was obtained 
from the Monterey Peninsula Airport observations during 1974. 
Hourly observations were examined to determine the signifi- 
cant mean daily wind vector for each day and corrections 
were applied to speed and direction when indicated by topo- 
graphic obstructions. Fortunately, no corrections were 
needed to northerly winds which are the driving mechanism 
for upwelling. 

An x-y coordinate system was established such that the 
surface wind stress Se = ole i. could be directly calcu- 
lated from the wind component parallel to the coast. The 
Ekman mass transport (S); surface current velocity ve 
and the average velocity within the Ekman layer (U) were 


calculated by the following relations: 


V2 29 ta ee (5b) 
O y 
Bo 





minere : V Surface current vector 


and Z 

—- iL 

U = 5 f va a 

O 

where: U = average velocity in x direction in the layer 
and: 

Sy = tL /E (6a) 
where S. = mass transport in x direction (offshore) 


x 


The vertical current representing upwelling is determined 
Per unit y (along coast) from the continuity equation, 


a” 


Cw vie= 0, such that 


UCD) = WC(L) 5) 
where: U = average horizontal velocity in the Ekman 
layer (m/day) 
W = upwelling speed (m/day) 
D = Ekman depth (m) 
L = horizontal width of the upwelling region (m) 
Zero current is assumed parallel to the coast. 

The horizontal width of the upwelling on (L) was 
determined from Yoshida's (1967) equation. Applying typical 
values to this equation yields a horizontal width of 3 X 
mer im. 

Principles of continuity were used to arrive at the mean 
daily upwelling current values and these figures were aver- 


aged over a 30-day period. 


34 





B. RADIATION INDEX 

peeewiscd ineident solar radiation index which takes in- 
Pemeocelm: the reduction of sunlight by local fog and stratus 
cloud cover was developed. Correlation of the measurements 
of surface irradiance and the strength of upwelling (as in- 
dicated by the upwelling index) by Tont (1975) is repre- 
sented by Figure 5. This figure is derived from the results 
ef Tont's study as shown in Figure 3 by plotting the percent 
of possible sunshine (Y-axis) as a function of the upwelling 
index (X-axis) which has been scaled to a maximum value of 
one. The resulting curves (A) and (B) show the conditions 
occurring after,and prior to, the upwelling maximum, respec- 
tively. - 

The difference exhibited between the curves may be due 
to changes in the character of the air mass brought about by 
seasonal migrations of the quasi-permanent thermal low over 
Nevada and the North Pacific sub-tropical high located west 
@& the coast. 

Figure 5 was used to calculate a revised incident radia- 
tion index by entering with the newly-developed upwelling 
index values and by applying the corresponding percent to 
the theoretical mean monthly radiation for clear sky condi- 
tions (possible sunlight). Figure 6 shows the possible sun- 
ieee throtighout the year «irom -ronc, 1975) as Qo? and the 
theoretical radiation at the surface (Q.) after the effects 


of fog and clouds are considered. 
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The incident solar radiation is an essential part of 
the ecosystem simulation model in that it is used in the 
equations governing the production of organic carbon by the 
photosynthetic activity of phytoplankton. These equations 


are discussed in detail by Pearson (1975). 


C. ALGAL SINKING er 

The phytoplankton sinking term which was added to the 
Pearson (1975) model is written in the form used by Riley 
(1965) in the units of the quantity of phytoplankton trans- 
ferred per day (g C7 See). The amount of phytoplankton 


that sinks out of the mixed layer ina day is given by: 
SNK = CV AZZ (26) 


where: SINK = flux of phytoplankton out of mixed layer 
(gC/m2-day) 


V = sinking rate (m/day) 


Z 


mixed layer depth (m) 
X2 = phytoplankton biomass in mixed layer (g/m?) 


V/Z = fraction of phytoplankton which sinks out of 
mixed layer per day 


Vertical circulation with the water column becomes 
Significant as the upwelling current speed and sinking rates 
of phytoplankton approach the same order of magnitude. Since 
the phytoplankton are non-mobile and are generally of the 
same density as the water or have developed shapes which 
increase their sinking resistance, they will be carried along 
With vertical currents in the water column. The vertical 


circulation during upwelling opposes sinking but may 
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accelerate downward transfer during brief periods of surface 
convergence and associated downwelling. 


The complete phytoplankton sinking term is: 
Sie = CCV — W/Z) XK? (27) 


where: W = upwelling speed (m/day) 

An average value of one meter per day was used’ for the sink- 
ing rate (V) in the computer simulation model. This value 
was determined from estimates presented by Lehman et al. 
(1975) and Bannister (1974). The SINK term acts to decrease 
the phytoplankton concentration in the modeled region and is 
subtracted in the phytoplankton equation (2) as shown: 


Grrr 
dt 





= PROD - RESP2 - GRAZ - SINK 2:3) 


2. Freer OrLANRTON ADVECTION “PERM 
An equation describing the horizontal advection of phyto- 
plankton is included in the revised model. This term is in 


Bmeminims @t flux and 1s written as: 


ADVEC2 W (X2)K G22) 


where: ° ADVEC2 = the change of concentration of phytoplank- 
ton over time (in the mixed layer) due to 
advection (g C/mé4-day) 


W = upwelling speed (m/day) 


X2 = phytoplankton concentration in the mixed 
layer (g C/m?@) 


=) 


K = advection coctavetent ann 
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The simplified advection term used in the model is de- 
rived from the equation (21) shown in the background theory 


eerion, which is: 


Advective change = W (X2) [- cease 


A Az! (21) 


oO 


The Ekman current offshore which results in an upwelling 
current (W) is effective over the depth of frictional resist- 
ance (D), but the phytoplankton in the model (which are 
uniformly distributed in the mixed layer) will be advected 
at the average velocity over the depth of the mixed layer as 
determined by the fraction of the mixed layer that is coinci- 
dent with the Ekman layer. It is known that the mixed layer 
varies seasonally from about 10 to 100 meters in depth and 
this causes an annual (seasonal) variation in the average 
horizontal current in the phytoplankton advection term. The 


average current of the mixed layer can be determined from: 


= , -7Z/D 


Z 
U = 7 f ¥,a2 where: Un =AePinermone of 
O 
G20) 

A shallow mixed layer experiences higher average veloci- 
ties than a deep layer as seen from the integral expression 
eq . “0D: 

A coefficient to describe the seasonal average of the 


vertically averaged horizontal current in the mixed layer is 


included in the advection equation as shown: 


Advective change = W (X2) [- a 


A m AZ! a 


8, 
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ire coefficient, Ls represents the seasonal average of the 
Fraction of the mixed layer which is coincident with the 
Ekman layer and is affected by the Ekman current. 

The constants in the advection expression (within the 
brackets of eq. 31.) are lumped into the advection coefficient, 
K, for convenience in equation (29). The range of K varies 


io Me me depending on the estimated values 


Erem about 10 
of the gradient distance, (AX) and (AZ) and the velocity co- 
efficient, Ka Aevalue oft 0:1 Ao was used in this model 


mor K. 


fee PREealLON LOSS 

The predation ae onee Function hypothesized for Monterey 
Bay is based on the assumption that the system supports a 
resident population of herbivorous zooplankton predators 
throughout the year. The pressure exerted on the herbivore 
prey by this eee level (N) of predators follows the 
curve in Figure 7 below the prey density threshold, d. The 
pressure is sufficiently low to allow growth of the zooplank- 
ton stock, but high enough to stabilize growth. A second 
pressure is imposed on the zooplankton (see Figure 7) after 
their density reaches the critical level, d. The simulation 
approximates conditions which may be brought about as tran- 
Silent predators move into the area presumably in response to 
increased food availability. In summary, predation pressure 
may be represented by a two part function with a pressure- 
density curve, a, due to resident predator species and curve 


b due to the transient predator population. 
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Figure 7. Predation pressure on zooplankton by 


transient predators. 
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Predator biomass may be expressed as shown in Figure 8, 
as a simple step increase or some other function of prey 
density. Consideration of the step type function 1s impor- 
tant because it allows the predator population (and conse- 
quent pressure) to maximize during those seasons when the 
zooplankton "standing crop" will support additional numbers 
of predators; the high predator pressure will deplete the 
Food stock to a lower level than would be reached by the 
resident predator pressure alone. Measurements of zooplank- 
ton biomass for 1974 (Traganza, 1975) suggest a rapid decline 
in "standing crop" following an early summer maximum. One 
might suspect that a transient predator population is forced 
out of the region once the food sources are depleted. This 
is done in the model by decreasing the predator population 
when the zooplankton biomass declines to a specified level. 
The predator lévels were related directly to zooplankton 
"Standing erop"” by triggering an increase in predator popula- 
tion at a threshold of zooplankton biomass during periods 
when this prey population was on the upswing. The predators 
were Similarly reduced at a second threshold during the de- 
clining phase of zooplankton growth. A set of four condi- 
flonase atements is used in the simulation to describe the 
predator-prey relationship. These statements specify the 
=O LitewaMic : 


tr 243 5 g AND x3 < 1.0 THEN: FISH = 1.0 


i 
awn 
© 


Peeks > 0 8 6ANDE see cee 
ae dt 
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TF == < 0 AND X3 > 0.2 GHEN: FISH = 5.0 
dX3 
i oe <5)  ANer 0.2 THEN: FISH = 1.0 


The modified predation pressure term in the zooplankton 


eauation of the Model is: 


HOSS = "Gi (F ISH) (X32) 37) 
where: 

LOSS = the amount of herbivorous zooplankton biomass 
lost per day as a result of predation 
(g C/mé day) 

L = loss rate or percent loss per FISH per day 

(day~1) 

FISH = number of predators (non-dimensional) 


Mo = goopitenkton "sitrandingwenep™ (Ce C/m*) 
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IVa REoULaES 


A. UPWELLING CALCULATIONS 

| The mesoscale wind data used to calculate upwelling for 
Eris Samiulation a shown by Figure 9. The plot depicts the 
frequency of occurrence (f = N/IN) for 45-degree segments 
of the compass. The prominent wind direction is as usual 
from the northwest .during the April to September upwelling 
period and shifts to the south in the first and last quarters, 
but this year there was a significantly high frequency of 
northwesterly winds during the last months of the year. The 
computer generated upwelling index (Figure 10) shows a pri- 
mary upwelling maximum in May, a second peak in about mid 
September and another in mid November. The possibility of 
a secondary upwelling peak is suggested by temperature obser- 
vations by Traganza et al. (1976). Figure 11 shows a rise 
in the isotherms peaking 30 to 45 days after the indicated 
wind initiated upwelling maximums in May and September. Al- 
though the water column does not respond instantaneously to 
surface wind stress, the delay noted here is excessive 
(Barham, 1957). The delay may be attributable to the assump- 
tion that the winds measured at Monterey Airport are actually 
representative of the winds over the bay, when in fact they 
may not be. The delay may also be caused by the lack of suf- 
ficient data to accurately depict the seasonal trends of 


temperature in the water column. 
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Traganza et al. (1976) observed a sharp rise in phosphate 
(Figure 12) with the major upwelling in May and a slight in- 
crease in phosphate and salinity (Figure 13) in September. 

No salinity data were available for the first half of the 
year. There was no clear-cut correlation of either phosphate 
or salinity to the upwelling index during November. Nutrient 
data from four to nine stations taken during seven cruises 

in 1974 were averaged (see Figure 12) and show a rapid rise 
ma OCtober which is probably not related to upwelling... "From 
the five year study of Monterey Bay by Bolin (1964) it is 
Keown that the nutrient concentration is characteristiecally 
low over the depth of the euphotic zone (0-200 meters) for 
two or three months at the end of the year. The restoring 

Oe thervnutrrente level Gia ail but the upper 20 to 30 meters 

to half of its May upwelling value while the surface tempera- 
ture reached a peak, may have signalled the beginning of the 
Davidson current period (see Bolin and Abbott, 1962). The 
Davidson current brings a southerly winter oceanic water 

mass into the Monterey region. This "Davidson water" is 
characterized by lower surface temperatures and surface salin- 
ities (due to high amounts of rain). The deeper water of 

the euphotic zone has higher salinities and nutrient concen- 
trations than exist at Monterey. The mechanism which estab- 
lishes these characteristics may be winter storm mixing oc- 
curring south of Monterey or upwelling and mixing from below 
brought about by divergent (cyclonic) eddies formed in the 
current stream as it moves northward along the coast 
(Geeeanzavet al., 1976). 
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Pagure 13. Seasonal salinity vartatron off Momeerey 
Bay. Values are the mean of four to nine 
stations (Traganza et al., 1976). 
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Traganza's (1976) data were obtained at stations in an 
X pattern referenced to a drogue (Figure 14) and at inter- 
vals of about two miles, two and one-half hours apart. The 
data from four to nine stations were averaged to obtain the 
values shown in Figures 11, 12 and 13. 

Barham (1957) argues that the vertical circulation of 
Monterey Bay during the upwelling period is characterized 
by a region of surface divergence coincident with the head 
of the Monterey submarine canyon. The effect is due to 
"channeling" of upwelled replacement water along the canyon 
axis with a plume appearing at the surface near shore (see 
Peeure 15). Evidence in support of this circulation concept 
is given by a distinct gradient of phytoplankton concentra- 
tion in the surface waters outlying the canyon head (Barham, 
1957). Barham believes that the gradient is due to the 
rapid horizontal advection of phytoplankton innoculum repre- 
senting a potential bloom away from the canyon head before 
the population shows any significant growth. During periods 
or reduced upwelling as indicated by the trends of the iso- 
therms, phytoplankton counts over the canyon rose to higher 
values, indicating that surface divergence or horizontal ad- 
vection had diminished. This possibility should not alter 
the timing of the model which is related to the general up- 
welling, but it does suggest serious spatial sampling prob- 


lems. 


Dg 








Figure 14, Chart of the Monterey area showing station 
location. 
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be VemilslLTEVITY ANALYSIS 

The model was tested under various conditions of phyto- 
plankton sinking rates, predation pressure on herbivorous 
zooplankton, and advection of phytoplankton to determine the 
effect on the simulation of the seasonal cycle phosphate, 
phytoplankton and herbivorous zooplankton. 

Losses of phytoplankton by sinking and therefore food 
limitation on zooplankton was examined by setting the preda- 
tion terms on zooplankton to zero and varying the rate of 
Baiemas OF algal cells, i.e. the availability of food. Sink- 
ing rates of zero, three, and six meters per day were used 
(Parsons and Takahashi, 1973 and Riley, 1965). The effects 
eee shown in Figures 16, 17, and 18. Under conditions of 
zero sinking, a rapid rise in zooplankton "standing crop" 
in late summer reduces the phytoplankton biomass quickly to 
a minimum on day 215. The large decrease in phytoplankton 
allows the nutrient concentration to remain at a relatively 
high level (Figure 16) since the uptake of nutrients by 
phytoplankton is decreased while ongoing zooplankton excre- 
tion and mixing by upwelling add to the nutrient concentra- 
tion. The rate of growth of the herbivorous zooplankton 
Popm@lataon or "standing crop" is shown to decrease with pues 
gressively higher phytoplankton sinking rates (Figure 18) 
as does the magnitude of the maximum zooplankton biomass 
reached during the year. The maximum also occurs later with 


increased sinking rate. 
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The sinking rate was varied in a second test but preda- 
Peep eessure was exerted on the zooplankton by employing 
the step function predator-prey relationship detailed earli- 
er. Additional predators were introduced in the simulation 
at a zooplankton biomass level of 1.0 g C/m* and were 
removed when zooplankton biomass fell to 0.2 g C/mé . Thresh- 
old values were determined empirically. The zooplankton bio- 


{x6 Ol Wet e C/m- 


mass maximum is reduced from 1-10 g C/m 
by the addition of predation pressure (Figure 21). The 
phospnate curves (Figure 19) show a marked response to phyto- 
plankton growth as evidenced by a lower nutrient minimum 
Following the peak of the summer bloom. A considerable 
change in the zooplankton response, 1.e. lower and earlier, 
1s evidently due to reduction of food sources as phytoplank- 
Gon 1S allowed to sink out of the mixed layer. The initial 
effect of increasing the sinking rate from zero to three 
m/day on zooplankton is a shift in the occurrence of the peak 
to approximately 70 days later. Further increase of the 
Sinking rate to six m/day results in a twofold decrease in 
the carbon biomass of zooplankton with an additional 20 day 
delay of the maximum. It 1s apparent that the rate of growth 
of zooplankton is slowed down and the maximum biomass occurs 
tater because of food limitation but that predator pressure 
is responsible for limiting the magnitude of the zooplankton 
biomass. The zooplankton peaks also occur earlier when pre- 
dation is applied and also occur betome the pnytopiamceen 
Bewess (Figure 20). The simulation, ~hemetore requiresea 
positive sinking rate. 
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The advection rate of phytoplankton was next varied 
while the sinking rate was held constant at 1.0 m/day 
(Bannister, 1974) and predator pressure was imposed in the 
same manner (and with the same thresholds) as in the second 
Mest. The results are shown in Figures 22, 23, and 24. A 
Slower rate of growth of zooplankton is produced again when 
mme advection coefficient, K, is increased from 0.0, to 0.3 
-_. Increasing the value of K has the direct effect of in- 
creasing the rate of phytoplankton advection. The single 
maximum of zooplankton simulated under conditions of high 
advection where K = 0.3 m > is due apparently to the sensi- 
“Gevity of the model's zooplankton growth equation to the 
mixed layer temperature and which permits a rapid growth 
of zooplankton during the mixed layer temperature minimum. 
Tne single peak of zooplankton (Figure 24) coincided with a 
me lacive Gemperature minimum about day 260. in Peamson’s 
Simulation a ten percent decrease in mixed layer temperature 
had a marked effect on zooplankton growth. This effect 
carried over to the partially modified simulation model. 

The effect of rapid zooplankton growth in a low temperature 
mixed layer, despite contradicting factors of growth such 

as low food availability, is due to the sensitivity of the 
respiration term in the zooplankton growth equation to tem- 
perature. This condition is a peculiarity of the model and 
it is doubtful whether the real ecosystem behaves in the 

same manner. The zooplankton peak precedes the phytoplankton 
peak when advection is zero indicating the model needs an 


g@avection term. 
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The predation pressure terms in the zooplankton equation 
mere varied (Figures 25, 26, and 27) while keeping the sink- 


ing rate at 1.0 m/day and the advection coefficient at K = 
ik 


O-lm~. The conditions of the test were: 
1. predation pressure set to zero; 
2. predation pressure defined by the linear function of 
zooplankton biomass: predator pressure = L times 
zooplankton biomass, where L = Q.01, and the units 


ae peC/m--day = (days !)(e C/mou. 
3.) wwe linea function in 2, but L = a7. 


4. the step function discussed under the methods section 
of this thesis. 


In all cases, the resulting zooplankton growth appears food- 
limited until about day 120 when phytoplankton growth begins 
to show a rate change. Complete removal of predator pressure 
allows the zooplankton to increase rapidly until food limita- 
tion (be depletion of phytoplankton) again occurs on day 275. 
When the predator pressure increases linearly (with L = 0.01), 


a Single zooplankton peak of 3.05 g C/m* e@eurs about day 


300. This is not consistent with observed conditions 
(Traganza, 1976). Doubling the predation pressure rate to 
Ea- © .02 day™* severely restricts zooplankton growth but pro- 


duces a curve with a hint of temporal conformity to the 
observed zooplankton maximum and minimum biomass levels. 
When the predator function is triggered at thresholds of 
zooplankton biomass, the zooplankton exhibit a double peak 


which is suggested by Traganza's (1976) data. 
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Socotra SON OF SIMULATION RESULTS WITH OBSERVED DATA 

The computer simulation was run with the new upwelling 
and radiation indices and the step function predation pres- 
sere term detailed in the preceding section. The sinking 
rate was set at 1.6 m/day and the advection coefficient set 
to K = 0.1m". 

A comparison of the simulation results with data observed 
me imaganza et al. (1976) is Shown in Pieures 26. 22 j aude 7c. 
The calculated nutrient values lag the averaged mixed layer 
concentrations by 20 to 30 days during the summer months with 
a moderate amount of error in magnitude. The nutrient maxi- 
mum occurs on day 144 with a value of 2.08 wg-at P/1 in the 
simulation. A rapid decrease in phosphate levels coincident 
With the summer phytoplankton bloom displays the impact of 
nutrient uptake by phytoplankton in the model. 

The simulated phytoplankton peak occurs 40 days after 
the nutrient peak during the period of maximum incident radia- 
tion (see Appendix C). The modeled response of phytoplankton 
was relatively inflexible with regard to the radiation index, 
me. the phytoplankton peak was follmd to GccUr within a Few 
days of the time of maximum incident radiation when the sink- 
ing rate was varied from about 1.0 m/day to 1.8 m/day. The 
relationship of the phytoplankton peak to the radiation peak 
is explained by examination of the phytoplankton growth equa- 
tions (shown in Appendix B). A slight increase in phytoplank- 
ton biomass was simulated during December (day 330) which 
could be due to the increased phosphate concentration at the 


end of the year. 
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Twenty days after the phytoplankton maximum, the model 
Simulated a zooplankton peak with a concentration of 0.97 
g C/m* . The zooplankton biomass growth corresponds well 
with the observed data until approximately the beginning of 
September. The measured values increase to a maximum of 
inmoo & C/m? in mid-December while the model simulated a 


recovery peak of 0.47 g c/n and slightly seaitver. 
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V. DISCUSSION 


The Monterey upwelling ecosystem simulation developed in 
this thesis has been shown to follow the seasonal trends of 
the observed phosphate, phytoplankton, and zooplankton data. 
There are errors in magnitude of the simulated response at 
various times of the year. The differences between observed 
and simulated zooplankton are most likely due to the fact 
that the model simulates herbivorous zooplankton while the 
observations include both herbivores and carnivores. Quanti- 
tative data defining the seasonal ratio of herbivorous to 
carnivorous zooplankton are needed to verify the simulation 
results. The encouraging aspect of the model is its ability 
to follow the seasonal trends and the fact that the timing 
of the response of phytoplankton biomass to the nutrient con- 
centration and the response of zooplankton to the phytoplank- 
ton cycle is biologically sound, 1.e. a moderate delay is 
Simulated between the peaks of successive trophic levels. 

Since all the factors likely to affect the three state 
variables (phosphate, phytoplankton and zooplankton) are not 
included in the model, some error in magnitude should appear. 
In recalling the philosophy expressed earlier, the objective 
of this thesis has been to create a time simulation of the 
dynamics of phosphate, phytoplankton and herbivorous zoo- 


plankton in a limited area over a long time. This objective 
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has been met by balancing the four neae characteristics of 
generality, realism, precision and simplicity. 

From studies by Barham (1957) and Bolin (1964) it is 
noted that there may be characteristic seasonal patterns of 
plankton dynamics in response to characteristic patterns 
of the hydrography of the Monterey upwelling ecosystem. It 
1s perhaps important to realize that the apparent limita- 
tions of the model as developed this far are relatively 
small when the use of the simulation to reproduce the 
ieeemeteri Stic patterns..1S..consa.deredypemuse:. 

As more and better data with which to define the dynamics 
are obtained, the simulation model of the Monterey upwelling 


ecosystem can progress to a more refined level. 
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VI. SWECESHT ONSGSPOReT URI HER RESEARCH 


The following comments are made as a result of questions 
which arose during the modification of the Monterey ecosystem 
model. This list is provided to identify some areas which 
may warrant further investigation. It is by no means com- 
plete, but should serve as a guide. 

1. The most significant outcome of this research has 
been the recognition of the importance of accurate forcing 
functions. The presence of the submarine canyon in Monterey 
Bay and the likely effect on upwelling suggests further 
effort be directed to verifying the upwelling index. 

2. The Monterey Bay area 1s not homogeneous in physical, 
chemical and biological parameters as evidenced by the 
patchiness of biological samples and the observable spatial 
gradients of the physical and chemical properties. 

3. Additional work on the biological aspects of the 
model including verification and refinement of the predation 
terms and development of an appropriate kinetic expression 
Pome natea ted . 

u, The simulation appears overly-sensitive to tempera- 
ture variations in the mixed layer. 

5. Various investigators have shown the importance of 
phytophagous fish in ecosystem models. A fish herbivore 


term should be investigated. 


Td 





6, Sole question as to the stability of the CSMP 360 
routine over long time simulation has been raised. This 
might be investigated and resolved by IiNdetia sealer ene 
model part way through the year. 

7. The effect OF varying the forcing functions in the 
revised model has not yet been studied. 

8. The advection expression might be expanded to in- 
clude the effects of all current systems influencing the 


Monterey region. 
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